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Nederlands 


Deze rapportage is een verslag van onderzoek dat is uitgevoerd in het kader van fase 1 van 
een onderzoeksproject door het CBS in opdracht van Staatstoezicht op de Mijnen (SodM). 
Dit onderzoek is ten behoeve van een statistische onderbouwing van het meet- en 
regelprotocol voor gasexploitatie in de provincie Groningen. Het onderwerp van dit rapport 
is een heranalyse van trends in de bodemdaling in de provincie Groningen, gerapporteerd in 
fase 0, in December 2014. Voor de voorliggende analyse is de tijdreeks voor de GPS 
gegevens aangevuld tot aan 3 Januari 2015. Daarnaast zijn alle GPS data opnieuw 
gecalibreerd, in een onafhankelijk onderzoek uitgevoerd door de TU Delft. 

Zoals ook bleek uit de eerdere analyses is er een statistisch significante trendbreuk in de 
bodemdaling ongeveer 2 maanden nadat de productie sterk was gereduceerd. De 
trendwijziging kan zich geleidelijk gemanifesteerd hebben over een periode van enkele 
weken, en er is daarom een onzekerheid van ruwweg een week of twee over de centrale 
datum van deze overgang. 

English 

This is a report on the research that has been carried within phase 1 of a project being 
carried out by Statistics Netherlands and commissioned by State Supervision of Mines 
(SodM). This research is part of the underpinning of the statistical methods employed to 
support the protocol for measurement and regulation of the production of natural gas in 
the province of Groningen. The subject of this report is to re-analyse the trends in the 
ground subsidence in and around that region, first reported in phase 0, in december 2014. 
For this re-analysis the time series of the GPS data has been updated with more recent 
measurements up to Jan. 3 2015. The entire GPS data set has also be re-calibrated, 
independently from the previous calibrations, by researchers at the TU Delft. 

In accordance with the earlier report, it is found that there has been a significant change in 
the ground subsidence. The changeover in the trend of ground subsidence can have 
become manifest gradually over a period of several weeks, which implies that the central 
date of the transition is also uncertain by roughiy a week ortwo. 


1 Introduction 


In this update for the most part the same methods are applied to the newtime series of GPS 
data. Forthis reasonthe reader is referred to the earlier report for an overview of the methods 
and the methodological background. The time series being analysed in this report contains the 
entire period that was analysed in the previous report, but is extended forward to coverthe 
epoch up to January 3 2015. Also, the most up-to-date Information is being used regarding 
reference stations. In practice this means that the entire set of GPS data is re-calibrated for all 
of the stations, and for the entire period. 
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Figure 1 The relative locations with names of each of the GPS stations from which 
data are availahle. The red diamond symhols are locations within the area of the 
production field 'Middle', where production has heen reduced. The stations indicated 
with hlue squares are locations outside of the main field, collectively designated as 
'Edge', although at the southwestern station Eemskanaal there are wells in 
production. The two stations indicated with green triangles are also production 
locations in a separate part of the field 'South' where production has not heen 
reduced. Well locations where production has heen reduced are shown as hlack 
crosses, with in addition three wells very close to this region (AMR, TJM, SDB). The 
indicated scales are in km 



2 The filtered time series for each station 


The first processing step of the time series of the GPS stations is to apply the filtering procedure 
with the filterfactors as derived in the previous report. As is done in the previous report, the 
average height of each station is subtracted. The results are shown in figure 2 

The quantity of interest is the 'sagging' of the stations in group 'Middle' and in group 'South' 
with respect to those in group 'Edge'. It is this differential displacement over a spatial scale of 
the order of, or smallerthan, the field that will lead to the build-up of stresses. Therefore the 
displacements, averaged for each group, are subtracted in the sense 'Middle' minus 'Edge' and 
'South' minus 'Edge'. The result is shown in figure 3. 

To assess the extentto which correlations might still be present in the differential displacement, 
two additional differential time series are shown in figure 4. 

One measure concerns the differential displacement between stations that are all within group 
'Middle', in the sense (Zeerijp+ten Post) - (Stedum+Overschild) which is shown as an orange 
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Figure 2 The time series of the GPS height, after filtering out intermediate time scale 
variations. For presentational purposes for each station the average of its time series 
over the period from day 50 to day 214 is suhtracted. The vertical scale is in m. Top 
Panel: GPS stations near the centre of the production field, hottom panel: GPS 
locations in group edge. 
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line. The other concerns the differential displacement between stations that are all within group 
'Edge', in the sense (Usquert+Veendam) - (Delfzijl+Eemskanaal) shown as a purple line. In the 
absence of remaining biases these differential measurements should not show any trends. 
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Figure 3 The time series of the filtered GPS height for the differences hetween the 
group averages 'Middle' — 'Edge' (in red), and 'South' — 'Edge' (in yellow). 



While it is clear that for most of 2014 no trend is present in these sets of differential 

Figure 4 The differential time series of the filtered GPS height within groups. Purple 
line: differential measurements hetween stations all within group 'Edge'. Orange line: 
differential measurements hetween stations all within group 'Middle'. 
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displacements, in Octoberthere is a sudden change in particular in the group 'Edge', but also to 
some extent in the group 'Middle'. At this stage it is unclear whetherthis is caused by further 
calibration issues that require resolving, orwhetherthe altered spatial pattern of production 
could have influenced these measures. 
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3 Drift trends 


Since the reduction of the gas production at some locations near stations in the group 'Middle' 
has been rapid, it is of interest to explore whether a break can be found in the linear component 
of the downward drift forthe time series for'Middle' minus 'Edge'. To this end one can fit not 
only a single straight line, using Standard least-squares fitting, but also introducé a break-point 
with a different straight line fit before and afterthat point. Given thatthe latterfit has more 
degrees of freedom than the former, this should be accounted for in assessing whetherthe 
two-section fit with a break is genuinely betterthan a single straight line. The Akaike 
Information Criterion (AlC) is an equivalent measure to the reduced measure used in the 
previous reports. 

Figure 5 AIC score as a function of the position of a break in the linear model fit of 
the difference between the group averages of'Middle' and 'Edge'. The minimum for 
'jump' is at 2014-02-06 and that of'kink' at 2014-04-02. 
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Figure 6 Filtered times series for the group averages of 'Middle' and 'Edge' with the 
fitted trends. 



4 Alternative model 


One possible issue with the previous fit of the break in the trend, is that the number of 
measurement contributing to each of the average lines changes with time. Before Februari 
2014 there is only one measurement station contributing to each of the lines. The previous 
model does not take this into account: one would expect the model to give more weightto 
differences at later dates as more measurement stations are contributing to the averages. 

One possibility to avoid this problem is to not calculate the difference, hut to model all time 
series combined. The difference in slopes between the two groups can be modelled explicitly 
and it is possible to test if this difference is significant. 

A piecewise linear continuous fit can be written as 

Vi = Po + Pik + Pisit - to^ti + ei, (1) 

with s(t) the step function {s(t) = 1 for t > 0 and 0 otherwise). We can then fit each of the time 
series with a piecewise linear continuous fit: 

^{Po+Pik +P 2 s(t-to)ti + ei if'Edge', 

{(.Po + Ps) + (.Pi + P4)ti + (.P2 + P5)^(.t - to)^i + if'Middle'. 

We can then test if the coëfficiënt P^ is zero. In figure 7 the coëfficiënt significantly differs from 
zero. Also in 'Edge'the discontinuity is significant {beta^ A 0). Therefore, both 'Edge' and 
'Middle' decrease in time and for both this decrease decreases from approximately 2014-03-14. 
However, for'Middle'this decrease is much strongerthan for'Edge', suggesting that this 
difference is related to the gas production. 
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Figure 7 Filtered time series of each of the measurement stations in 'Edge' and 
'Middle' with the fitted model. The hreak is located at 2014-03-14. 



5 Re-processing of GPS data by TU Delft 


The GPS data used has been processed from the original satellite measurements using one 
Standard algorithm and associated method. In orderto assess some of the assumptions 
regarding sources of errorand their effect, a different algorithm has been applied, developed at 
the TU Delft. Their dataset of reduced GPS signal per GPS ground location has been made 
available, so that the same difference technique can be applied, as has been done forthe 
previous reduced GPS datasets. The data are shown in figure 8, which is the equivalent of figure 
2 . 

Fig. 9 shows the 'sagging' of the stations in group 'Middle' and in group 'South' with respect to 
those in group 'Edge', ie. the equivalent of what is shown in figure 3. 

The differential displacement between stations that are all within group 'Middle', in the sense 
(Zeerijp-Ften Post) - (Stedum-FOverschild) is shown as an orange line in figure 10. The other 
measure concerns the differential displacement between stations that are all within group 
'Edge', in the sense (Usquert-FVeendam) - (Delfzijl-FEemskanaal) shown as a purple line. In the 
absence of remaining biases these differential measurements should not show any trends. 
Comparing figure 10 with figure 4 shows that forthe differential time series between the 
stations that are all in the group 'Edge'the behaviour after October 2014 is very different. This 
would appearto indicate that in that epoch there are indeed some calibration issues in the 
earlier dataset which the reprocessing by the TU Delft have resolved. The same epoch does still 
show some anomalies forthe time series involving stations within the group 'Middle'. 
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Figure 8 The time series of the GPS height as determined hy the TU Delft method, 
and after filtering out intermediate time scale variations. Note that the vertical scale 
is in mm. Top Panel: GPS stations near the centre of the production field, hottom 
panel: GPS locations in group edge. 
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5.1 fitting to the TU Delft data without additional filtering 

While the fitting procedures used are robust to the precise properties of the errors, it is 
worthwhile to explore to what extentthe point-to-point correlation introduced by the filtering 
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Figure 9 The time series of the TU Delft processed GPS height data, and after 
filtering, for the differences hetween the group averages 'Middle' — 'Edge' (in red), 
and 'South' — 'Edge' (in yellow). Note that the vertical scale is in mm. 



Figure 10 The differential time series of the TU Delft processed GPS height data, and 
after filtering, within groups. Purple line: differential measurements hetween 
stations all within group 'Edge'. Orange line: differential measurements hetween 
stations all within group 'Middle'. 
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process affects the conclusions drawn from the trend fitting. At this stage it is worth 
remembering that the choice of fits does not have an underlying physical model forthe 
subsidence. The intention is merely to assess whetherone can rejectthe hypothesis that there 
is no change in average subsidence rate overthe entire period of observation. 

Figure 11 is the equivalent of figure 5 but now the fitting is done on GPS data as processed by 
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the TU Delft without any additional filtering by CBS. While there are differences in the AlC for 

Figure 11 AIC score as a function of the position of a break in the linear model fit of 
the difference between the group averages of'Middle' and 'Edge'. The minimum for 
'jump' (red curve) is at 2014-04-05. For the 'kink' (yellow curve) there are several 
choices which are equally valid, which can be seen from the three broad minima. The 
final minimum on 2014-03-01 coincides within the uncertainties with what is 
obtained for the filtered data (fig. 5) 



Date 


the fits as determined forthe filtered and unfiltered data, as one would expect, in both cases 
one concludes that fits for which the average subsidence rate changes over the period of 
observation are superior and statistically significantly so. The data appearto pointto a trend 
change in late March early April 2014, with an uncertainty of a few weeks eitherway, possibly 
because the change is not quite as abrupt as the fits are. 


Statistics Netherlands | Scientific paper 2015 


12 











Figure 12 The fitted models, with either a jump or a kink, to the data as processed hy 
the TU Delft without any additional filtering for the hest-fitting position of that hreak 
in the linear model fit of the difference hetween the group averages of'Middle' and 
'Edge'. 



6 Conclusions 


Fromthe GPS data itcan be concludedthatthere is continued subsidence of the ground in the 
area of the wells where production was reduced in the month of January of 2014. Howeverthe 
rate of subsidence is lower after a date. The date of this break has been redetermined with the 
more extended and re-calibrated set of GPS data. In addition a completely independent data 
reduction of the GPS data, developed at the TU Delft, leads to the same conclusions in our 
analysis. The location of the break is around mid-March or very early in April, ie. about 9 weeks 
after the reduction in production, butthere is a fair margin of uncertainty (some 2 weeks) 
around the exact value of the time gap. 

While there is clear statistical evidence that a break has occurred in the subsidence rates, and 
the reduction in subsidence speeds is measurable with a high degree of precision, research is 
now in progress to perform a more systematic analysis of the correlation hetween the 
production time series and the GPS height time series. 
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